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Numerical  modeling  of  the  expansion  of  electric  thruster  plumes  provides  direct  means 
for  predicting  spacecraft  surface  contamination  and  erosion  due  to  plume  ions.  A  software 
package  named  COLISEUM  that  is  capable  of  self-consistently  modeling  plasma 
propagation  and  interactions  with  arbitrary  3-D  surfaces  is  being  developed  by  a  national 
team  of  researchers.  Despite  much  research  and  development  in  modeling  plume  expansion, 
it  is  necessary  to  continuously  validate  these  codes  using  laboratory  based  experimental  data. 
It  is  well-established  that  vacuum  chamber  facilities  affect  the  plume  of  these  devices.  Thus, 
the  models  must  not  only  describe  the  plume  expansion,  but  also  effects  of  the  vacuum 
chamber.  COLISEUM  has  been  designed  to  simulate  both  vacuum  chamber  configurations 
and  spacecraft  geometries.  This  work  provides  source  derivation  from  laser  induced 
florescence  (LIF)  data.  Included  is  a  study  that  compares  results  from  a  hybrid  particle-in¬ 
cell  model  (AQUILA)  with  Monte  Carlo  collisions  to  data  obtained  from  the  plume  of  Busek 
600 W  Hall  thruster  (BHT-HD-600).  This  data  includes  current  density,  velocity 
distribution,  and  energy  data. 


I.  Introduction 

Numerical  modeling  of  the  thruster  and  surrounding  environment  provides  direct  means  for  predicting  plume 
properties  where  experimental  methods  are  limited,  such  as  predictions  of  spacecraft  interactions.  Insight  into 
the  plume  properties  and  corresponding  spacecraft  interaction  would  provide  the  community  with  a  useful  tool.  The 
Air  Force  Research  Laboratory  (AFRL)  is  leading  the  development  of  COLISEUM  ,  a  3D  plasma  interaction 
framework  that  incorporates  a  plasma  expansion  tool  with  surface  interactions.  COLISEUM  has  been  designed  to 
be  usable,  flexible,  and  expandable.  COLISEUM  has  available  any  of  four  plasma  simulation  models,  RAY, 
prescribed _plume ,  DRACO,  and  AQUILA.  RAY  uses  a  ray-tracing  method  to  project  a  flux  from  a  point  surface. 
Prescribed _plume  imports  and  superimposes  a  plume  distribution  onto  surfaces.  DRACO  tracks  particles  along  a 
structured  Cartesian  mesh.  AQUILA,  the  focus  of  this  study,  is  a  hybrid  PIC  model  that  tracks  particles  along  an 
unstructured  tetrahedral  mesh.  COLISEUM  is  capable  of  modeling  chamber  effects,  which  are  known  to  affect  the 
plume  expansion2,  as  well  as  open  boundary  conditions. 
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Despite  much  research  and  development  in  modeling  plume  expansion,1  it  is  necessary  to  continuously 
validate  these  codes  using  laboratory  based  experimental  data.  This  paper  presents  a  study  that  compares  results 
from  an  AQUILA  simulation  to  experimental  data  from  the  BHT-HD-600  Hall  thruster.  The  numerical  study  uses  a 
simplified  Hall  thruster  to  simulate  an  electric  plume  in  a  chamber  environment.  Laser  induced  florescence  (LIF) 
data  is  used  to  construct  source  input  files.  Probes  in  AQUILA  collect  current  density,  velocity  distributions,  and 
ion  velocities  which  are  compared  to  experimental  data  for  code  validation.  The  full  capabilities  of  AQUILA  are 
expressed  by  using  probes  to  collect  additional  information  about  the  plume. 


II.  AQUILA  code 

AQUILA,  a  hybrid  particle-in-cell  model,  has  been  developed  in  the  framework  of  COLISEUM.  AQUILA  uses 
an  unstructured  tetrahedral  mesh  to  define  surface  and  volume  geometries.  The  geometry  and  mesh  used  by 
COLISEUM  can  be  produced  using  available  commercial  modeling  and  meshing  packages.  The  mesh  can  be 
loaded  into  COLISEUM  using  any  number  of  standard  forms  including  ANSYS  and  ABAQUS.  Individual  surfaces 
are  specified  to  distinguish  surface  properties  and  define  particle/surface  interactions. 

In  most  places  of  the  plume,  the  plasma  is  defined  to  be  quasi-neutral.  Following  this  assumption,  the  potential 
can  be  calculated  by  using  the  inverted  Boltzmann  equation: 


.  ,  kT 

=  +— In 


f  \ 
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\neoJ 


where  neo  is  a  reference  ion  density,  (j)o  is  a  reference  electrostatic  potential,  and  Te  is  the  reference  electron 

temperature.  For  plumes  where  the  quasi-neutral  assumption  does  not  hold,  such  as  behind  a  plume  shield, 
AQUILA  contains  a  non-neutral  solver. 


Collisions  are  performed  in  AQUILA  using  a  Direct  Simulation  Monte  Carlo  method  .  Collision  cross  sections 
for  elastic  collisions  between  neutrals  and  ions  are  calculated  using: 


cr 


el 

Xe-Xe 


2.1 17xl0“18 


v 


0.24 

rel 


el  _  el 

(JXe-Xe+  ~  ° Xe-Xe++ 


8.2807x10  16 


where  vrei  is  the  relative  velocity  between  the  two  particles.  Charge  exchange  collision  cross  sections  between 
ions  and  neutrals  are  defined  as: 


=  2.241 5x1 0“18  -2.7661xl0-19log(vre/) 
axe-xe"  =  l*10-2°(35.006  -  2.70381og(vre/))2 


The  current  density  probe  in  AQUILA  samples  particles  on  a  hemisphere  at  a  user  defined  location.  The  results 
along  the  hemisphere  are  averaged  to  give  values  from  0  to  90°. 

To  decrease  computational  time,  acceleration  techniques  have  been  introduced  into  AQUILA.  One  of  the 
acceleration  techniques  utilizes  a  subcycle  routine  that  decouples  the  ion  and  neutral  movements.  This  method  has 
been  demonstrated  previously  .  Subcycling  cycles  the  fast  particles  for  a  specified  time  step  before  moving  and 
injecting  the  slow  particles  and  performing  collisions.  Subcycling  in  AQUILA  is  shown  to  decrease  the 
computational  time  for  the  simulation  to  reach  steady  state  by  77%,  as  shown  by  Marshall  and  VanGilder  . 
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However,  it  is  also  shown  that  too  many  subcycles  (-100)  allow  an  ion  to  leave  the  main  region  of  the  plume 
without  undergoing  any  collisions,  an  event  that  may  not  be  entirely  physical. 

III.  Source  Definition 

Plasma  modeling  within  COLISEUM  begins  with  source  definition.  Sources  are  used  to  introduce  particles  into 
the  simulation  domain.  COLISEUM  provides  several  different  source  models,  all  specifying  a  velocity  distribution 
function  (VDF)  and  a  mass  flow  rate.  While  several  source  models  are  available,  the  test  case  used  in  this  paper 
used  the  source  FLUXRVZVR  .  As  the  name  suggests,  the  FLUXRVZVR  model  represents  the  exit  plane 
of  a  Hall  thruster  in  terms  of: 

1 .  particle  flux  versus  radial  position,  r 

2.  axial  velocity,  vz,  versus  radial  position,  r 

3.  radial  velocity,  vr,  versus  radial  position,  r 

Before  starting  the  simulation  run,  the  flux  function  is  converted  into  a  cumulative  distribution  function  for  radial 
position.  During  injection  of  particles,  the  simulation  uses  the  thruster’s  mass  flow  rate  to  determine  the  number  of 
particles  to  create  at  each  time  step.  The  previously-computed  cumulative  distribution  function  (CDF)  is  then  used  to 
place  source  particles  at  radial  distances  with  the  correct  probability.  Once  the  radial  location  of  the  injected  particle 
is  known,  the  initial  velocity  components  are  calculated  from  the  velocity  functions.  In  addition,  thermal 
components  are  added  to  these  velocities  based  on  user-specified  temperatures. 

Laser-induced  florescence  data  provides  a  natural  way  of  determining  the  VDF  required  by  the  simulation  source 
model.  After  processing,  the  data  taken  at  each  location  gives  a  probability  distribution  function  of  the  velocity 
component  aligned  with  the  laser  used  to  probe  the  plasma.  To  get  velocity  distribution  functions  along  a  different 
axis,  the  orientation  of  the  laser  is  changed.  For  this  source  model,  all  VDF  inputs  were  taken  from  LIF  data  at  an 
axial  distance  of  15  mm  downstream  from  the  thruster  exit  plane.  Figure  1  shows  the  geometry  of  the  BHT-HD- 
600. 


Figure  1:  BHT-HD-600  geometry  for  LIF  analysis.  The  data  for  the  LIF  analysis  is  taken  at  y  =  15  mm  which 

is  between  the  exit  plane  of  the  thruster  and  the  cathode. 
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A.  Velocity  Distributions 

Figure  2  shows  a  representative  axial  VDF  for  the  BHT-HD-600  taken  at  a  radial  distance  of  24  mm  from  the 
centerline.  A  distinct  velocity  peak  is  seen  around  20,000  m/s.  To  quantify  the  mean  and  spread  of  the  peak,  the 
velocity  distribution  function  is  converted  to  a  histogram  and  a  Matlab  function  that  fits  a  specified  number  of 
Gaussians  is  used.  Figure  3  illustrates  this  process  -  the  red  lines  show  the  fitted  Gaussians.  In  most  cases,  two 
Gaussians  were  used  to  fit  the  histogram  -  one  captures  the  property  of  the  peak  of  interest  while  the  other 
represents  background  noise  in  the  signal. 


Figure  2:  Axial  velocity  distribution  function  at  r  =  24mm  of 
BHT-HD-600. 


Figure  3:  Histogram  with  fitted  Gaussians  for  axial  velocity 
distribution  function  at  r  =  24  mm. 
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A  similar  procedure  is  followed  at  each  radial  location  with  LIF  data  for  both  the  axial  and  radial  velocity 
components.  The  mean  of  the  main  Gaussian  is  used  to  represent  the  velocity  magnitude,  while  the  standard 
deviation  is  used  to  represent  the  temperature.  Since  the  source  model  assures  azimuthally  symmetry,  only  values 
from  one  side  of  the  thruster  are  needed  -  for  this  case,  values  from  the  non-cathode  side  of  the  thruster  are  taken 

Upon  further  inspection  of  the  LIF  data,  it  is  observed  that  evidence  of  a  second  ion  population  is  seen  in  the 
axial  velocity  distribution  functions  found  on  the  cathode  side  of  the  thruster.  Compared  to  the  above  figures,  Figure 
4  shows  the  corresponding  axial  distribution  at  r  =  24  mm  on  the  cathode  side. 


Figure  4:  Axial  velocity  distribution  function  at  r  =  -24mm  (cathode  side)  of  BHT-HD-600. 

Two  defined  axial  velocity  peaks  are  observed,  but  when  the  radial  distributions  are  investigated,  no  difference 
in  behavior  is  seen.  It  is  conjectured  that  the  second  axial  peak  on  the  cathode  side  may  be  related  to  additional 
ionization  that  occurs  near  the  exit  plane.  Electrons  streaming  from  the  cathode  are  assumed  uniformly  azimuthally 
distributed  by  the  magnetic  field  once  they  enter  the  thruster  acceleration  channel  and  as  a  result,  ionization  is  also 
uniform.  However,  neutrals  outside  the  thruster  on  the  cathode  side  are  more  susceptible  to  ionization  by  electrons 
before  they  reach  the  channel.  Since  these  ions  are  formed  outside  of  the  main  accelerating  potential,  their  existence 
is  manifested  by  a  smaller  lower-energy  peak  in  the  axial  distributions.  As  currently  written,  the  AQUILA  source 
model  cannot  model  an  asymmetric  exit  plane  distribution.  Nevertheless,  source  model  information  including  this 
second  population  is  also  processed  in  case  one  considers  the  second  low-energy  peak  important. 

B.  Flux  Distribution 

Flux  measurements  for  the  BHT-HD-600  have  not  yet  been  completed.  An  attempt  was  made  to  extract  ion  flux 
data  from  LIF  signal  strength.  Signal  strength  in  the  linear  regime  is  proportional  to  both  ion  density  at  the  probed 
state  and  to  the  intensity  of  the  beam  .  This  approximation  holds  if  the  ion  and  electron  temperatures  and  the 
electron  density  are  relatively  constant  throughout  the  region.  The  peak  signal  strength  at  each  data  location  for  z  = 
15  mm  is  plotted  as  single  points  in  Figure  5.  Using  a  curve  fitting  algorithm,  a  Gaussian  profile  was  fit  to  this  data. 
The  flux  profile  shifts  the  center  of  the  distribution  away  from  the  centerline  of  the  channel  (at  r  =  -28  mm)  towards 
the  center  line  of  the  thruster  (r  =  24  mm).  An  input  file  was  generated  using  this  distribution. 
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Figure  5:  Flux  profile  generated  from  LIF  signal  strength  data  fitted  with  a  Gaussian  distribution.  The 
Gaussian  fit  to  the  LIF  data  is  centered  at  r  =  -24  mm  rather  than  over  the  centerline  of  the  thruster  at  r  =  -28 

mm. 


IV.  Problem  Description 

The  chamber  and  thruster  geometry  for  the  simulations  is  based  on  the  BHT-HD-600  Hall  thruster  live  tests 
inside  Chamber  6  at  the  Air  Force  Research  Laboratory.  During  the  thruster  testing,  graphite  panels  were  added  to 
the  chamber  to  lower  sputter  rates  of  the  chamber  and  reduce  re-deposition  of  the  chamber  materials  back  onto  the 
thruster.  The  simulation  preserved  the  placement  and  orientation  of  the  graphite  panels  while  simplifying  the 
geometry  of  the  panels  and  other  components  to  increase  mesh  quality.  Figure  6  shows  the  surface  mesh  of  the 
chamber,  graphite  panels,  and  thruster  orientation.  The  thruster  fires  in  the  negative  z  direction. 

The  chamber  modeled  is  a  1.8  m  diameter,  2.9  m  length  stainless  steel  chamber.  The  graphite  panels  are 
modeled  as  semi-circular  shapes  concentric  with  the  chamber,  60  cm  in  width  and  lengths  from  79  cm  to  1.22  m. 
Cryogenic  pumps  are  also  included  in  the  geometry  at  the  rear  of  the  chamber,  87  cm  wide  and  83  cm  apart.  The 
chamber  background  pressure  was  set  to  7  x  10"4,  corrected  for  Xe. 
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Y 


Figure  6:  Chamber  6  geometry  and  mesh.  The  green  elements  are  the  graphite  panels.  The  blue  elements 
are  the  cryogenic  pumps.  The  thruster  is  the  orange  box  in  the  center  of  the  chamber. 


The  simulated  Hall  thruster  geometry,  shown  in  Figure  7,  is  the  plasma  emitting  source  in  the  simulation. 
Particles  are  emitted  from  the  annular  region  of  the  thruster  face.  Based  on  the  geometry  of  the  BHT-HD-600,  the 
annulus  has  an  inner  diameter  of  24  mm  and  an  outer  diameter  of  32  mm. 


Figure  7:  BHT-HD-600  simplified  mesh.  The  yellow  elements  are  defined  as  the  particle  emitting  source. 

Notice  the  change  in  element  resolution  from  the  front  face  to  the  back. 

The  volume  mesh  used  in  this  simulation  is  shown  in  Figure  8.  The  mesh  is  a  structured  mesh  in  the  core  region 
of  the  plume,  produced  to  increase/control  the  resolution  of  elements  in  the  immediate  region  of  the  plume.  This 
region  is  defined  by  the  annulus  extending  outward  30  cm  from  the  thruster  face,  with  the  end  elements  triple  the 
size  of  the  original  elements.  An  unstructured  mesh  fills  in  the  remainder  of  the  volume. 
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Figure  8:  The  simulated  volume  mesh  with  a  structured  mesh  in  the  region  of  the  plume  and  an 
unstructured  mesh  in  the  remaining  volume.  The  structured  mesh  is  used  to  control  the  resolution  of  the 

mesh  in  the  immediate  region  of  the  plume. 


The  annular  region,  shown  in  Figure  7,  is  the  particle  emitting  source  of  the  simulation.  Based  on  experimental 
findings15,  the  particles  emitted  include  not  only  Xe  neutrals  and  singly  charged  Xe  ions,  but  multiply  charged  Xe 
ions  as  well.  The  distribution  of  neutrals  is  defined  by  a  drifting  Maxwellian  source  model  with  a  thermal  drift 
velocity  of  297  m/s  and  a  temperature  of  700  K.  The  mass  flow  rate  of  neutral  particles  into  the  simulation  reflects 
10%  of  the  unionized  propellant.  Source  files  for  doubly  charged  ions  were  generated  using  the  singly  charged  ion 
source  files,  increasing  the  magnitude  of  velocity  by  the  V2  to  represent  an  increase  in  speed  related  to  their  charge 
according  to  the  energy  relation: 


—mv2 

2 


=  Zq(j) 


where  Z  is  the  ion  charge,  q  is  the  elementary  charge,  and  (/)  is  the  potential. 

Although  evidence  of  a  low  velocity  ion  population  was  seen  in  the  axial  LIF  data,  this  population  was  only 
evidenced  in  the  cathode  region  of  the  plume.  Since  the  majority  of  the  plume  lacks  this  second  population,  the 
simulation  was  completed  using  only  the  main  population  of  ions,  and  the  flux  profile  shown  in  Figure  5.  Ion 
temperatures  of  1 0  eV  in  the  radial,  axial,  and  azimuthal  directions  were  used. 

The  electrostatic  potential  uses  the  quasi-neutral  potential  solver  in  AQUILA.  The  reference  potential  (j)0  is  set  to 


20  eV,  and  the  reference  density  neo  is  sampled  from  the  domain  at  a  point  centered  on  the  annular  channel.  The  ion 
temperature  T  is  specified  using  a  polytropic  model: 


T  =  T 


f  V-1 

na 


\neoJ 


where  the  reference  electron  temperature  Te  is  set  to  5.0  eV  and  y  is  1.3. 


The  timing  scheme  in  AQUILA  is  used  to  determine  when  the  simulation  reaches  a  steady  state.  Subcycling  was 
used  in  this  simulation  to  decrease  computation  time.  The  number  of  subcycles  was  set  at  10,  with  an  ion  time  step 
of  2  x  10"7  s,  and  the  number  of  cycles  was  2500.  The  region  of  the  domain  specified  as  the  cryogenic  pumps  takes 
incoming  particles  out  of  the  domain  at  a  user  specified  rate.  To  determine  correct  cryogenic  pumping  speed  for  the 
simulation,  the  neutral  count  should  reach  a  steady  state  when  the  source  ions  begin  to  collide  with  the  wall.  Using  a 
time  step  of  2  x  10"7  s,  source  ions  reach  the  wall  in  7.5  x  10"5  s.  During  the  simulation,  this  occurs  after  39  cycles. 
Running  the  simulation  with  various  cryogenic  sticking  coefficients,  it  is  shown  that  this  is  achieved  at  a  pumping 
speed  of  8%.  Figure  9  shows  the  particle  count  plot  for  neutrals. 
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Cycles 

Figure  9:  Neutral  particle  count  for  different  sticking  coefficients.  An  adequate  sticking  coefficient  of  8% 
is  indicated  by  a  level  neutral  count  after  39  cycles,  when  source  ions  begin  to  collide  with  the  wall. 


V.  Comparison  against  Experimental  Data 

The  plume  studies  of  the  BHT-HD-600  conducted  by  Ekholm,  et  al  provided  measurements  of  the  ion  current 
density  profile,  ion  energy  distributions,  and  ion  species  fraction  distributions  using  a  nude  Faraday  probe,  retarding 
potential  analyzer  (RPA),  and  ExB  probe.  LIF  measurements  were  taken  from  Charles  and  Hargus  .  This  suite  of 
data  serves  as  a  comparison  to  the  test  cases  completed  using  AQUILA  and  the  input  parameters  defined  above. 

A.  Current  Density 

The  current  density  probe  in  AQUILA  was  set  up  to  sample  particles  along  a  hemisphere  60  cm  in  front  of  the 
thruster,  allocating  200  bins  for  sampling.  Figure  10  shows  the  current  density  from  the  test  case  compared  to 
experimental  results  from  Faraday  probe  measurements  at  60  cm.  As  shown  in  Figure  10,  the  trend  in  the  current 
density  profile  agrees  with  experimental  data  across  the  region,  although  the  magnitudes  only  agree  from  -20°  to  20°. 
In  the  wing  region  of  the  plume,  the  data  for  the  test  case  is  lower  than  what  is  seen  experimentally.  Experiments 
are  known  to  have  difficulties  in  this  region  due  to  charge  exchange  collisions  and  secondary  electron  emission  from 
ion  impact  with  the  probe.  The  difference  could  also  indicate  greater  collision  rates  in  the  plume  than  what  is 
modeled.  Increasing  collision  rates  would  lead  to  greater  plume  divergence  and  spread  more  ions  out  towards  the 
wings.  Despite  the  differences,  the  similar  trend  between  the  two  sets  of  data  suggests  that  only  minor  refinements 
are  needed  in  the  model  to  bring  the  both  sets  of  data  into  agreement. 
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Figure  10:  Current  density  at  60  cm  from  Ekholm,  et  al  and  from  test  case.  The  trends  in  both  sets  of 
data  are  similar,  with  the  greatest  differences  occurring  in  the  wings  of  the  plume  from  -40°  to  40°.  Greater 
collision  rates  in  the  model  would  broaden  the  spread  of  the  test  case  profile. 


B.  Near  Field  Ion  Velocities 

Figure  1 1  shows  a  comparison  of  axial  ion  velocities  at  x  =  0  mm,  z  =  60  mm  from  the  thruster  exit  with  LIF 
peak  data  values  taken  by  Charles  and  Hargus  and  an  uncertainty  of  500  m/s.  Figure  12  is  the  radial  velocity 
comparison  in  this  same  region.  AQUILA  tracks  particles  based  on  the  collisions  they  undergo  and  attaches  a  prefix 
to  the  species  name.  An  SRC  delimiter  signifies  a  beam  ion  that  has  not  experienced  any  collisions.  ELCOL  and 
CEXCOL  label  ions  that  have  undergone  elastic  and  charge  exchange  collisions,  respectively.  It  can  be  noted  that 
the  trend  of  the  simulated  source  ions  follows  the  trend  of  the  LIF  data,  with  the  axial  velocities  around  20  km/s  and 
a  3  km/s  deep  depression  around  y  =  0  mm.  The  radial  velocities  follow  a  linear  pattern  across  both  sections  of  the 
channel,  changing  from  a  divergent  to  convergent  profile  at  y  =  +-30  mm,  noted  by  a  shift  in  velocities  from  positive 
to  negative.  As  is  expected,  the  ions  undergoing  collisions  are  shown  to  have  axial  velocities  around  1000  m/s  and 
radial  velocities  near  0  m/s.  There  is  also  a  significant  amount  of  distribution  in  the  velocity  which  is  seen  in  data 
taken  by  Charles  and  Hargus.  While  the  greatest  population  of  ions  occurs  at  the  peak,  higher  velocity  and  lower 
velocity  populations  exist.  Figures  11  and  12  suggest  that  the  model  portrays  consistent  results  in  the  near- field  as 
compared  to  experimental  data. 
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at  x=  0  mm,  z  =  60  mm 


Figure  11:  Axial  LIF  comparison  at  60  mm  in  front  of  thruster.  The  trends  from  both  the  experimental 
and  test  case  data  agree,  with  a  significant  amount  of  distribution  in  the  test  case  data  also  seen  in  the  data  by 

Charles  and  Hargus. 


at  x=  0  mm,  z  =  60  mm 


Figure  12:  Radial  LIF  comparison  at  60  mm  in  front  of  the  thruster.  The  trends  from  both  the 
experimental  and  test  case  data  agree,  with  a  significant  amount  of  distribution  in  the  test  case  data  also  seen 

in  the  data  by  Charles  and  Hargus. 
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C.  Far  Field  Energy  Distribution 

The  energy  distribution  of  ions  in  the  simulation  can  be  determined  from  the  ion  velocities.  Using  the  equation 
for  kinetic  energy:  KE  =  Vi  mv2,  with  the  mass  of  a  Xe  ion  as  2.182  xlO"25  kg,  the  kinetic  energy  per  charge  is 
computed  along  an  arc  60  cm  in  front  of  the  thruster  in  the  x-z  plane.  The  resulting  plot  is  shown  in  Figure  11, 
compared  against  RPA  experimental  data.  The  RPA  probe  measures  an  energy  distribution  per  charge;  the  peak 
values  are  shown  in  Figure  11,  along  with  a  bar  representing  the  distribution  spread  seen  in  the  data,  around  10%  of 
the  peak  value.  The  same  particle  labels  are  used  in  Figure  13  as  in  Figure  1 1  and  12.  As  is  expected,  collision  ions 
are  shown  to  have  energies  between  0  and  50  eV/q,  with  some  elastic  scattering  seen  around  100  to  150  eV.  It  is 
interesting  to  note  that  Ekholm,  et  al  noticed  a  peak  in  the  RPA  data  in  the  150  eV  range,  attributing  it  to  elastic 
scattering.  Beam  singly  charged  ions  have  much  higher  energies,  with  a  mean  value  near  325  eV/q  and  a  broad 
distribution  from  150  eV/q  to  550  eV/q  that  falls  outside  the  range  of  experimental  data.  The  doubly  charged  beam 
ions  have  the  same  amount  of  energy  as  their  singly  charged  counterparts,  but  because  of  their  double  charge,  their 
energy  per  charge  is  only  half.  Notice  that  the  beam  ions  reside  only  in  the  region  between  -40°  to  40°,  suggesting  a 
low  beam  divergence.  The  combination  of  a  large  range  of  ion  energies  and  mean  value  higher  than  the  peak 
experimental  results  suggests  that  the  source  may  be  emitting  particles  at  velocities  that  are  too  high  and/or  over  too 
great  a  range.  By  adjusting  the  thermal  component  of  the  velocity,  it  is  possible  to  confine  the  velocity  range  and 
bring  the  model  energy  values  to  a  more  reasonable  level.  The  near  field  results  suggest  that  the  source  model 
velocities  are  consistent  with  experimental  FIF  data,  although  the  far  field  results  indicate  that  the  particle  velocities 
are  too  high.  This  suggests  that  there  may  be  an  issue  in  the  far  field  that  needs  to  be  addressed,  including  mesh 
resolution,  potential  solver  parameters,  collision  properties,  and  background  pressure. 


at  r  =  60  cm,  y  =  0  cm 


Degrees  off  Centerline 


Figure  13:  Ion  Energy  Distribution  at  r  =  60  cm  from  Ekholm,  et  al  and  from  the  test  case.  The  ions  are 
labeled  according  to  the  collisions  they  experience.  EL  COL  and  CEX  COL  denote  elastic  collisions  and 
charge  exchange  collisions  respectively.  SRC  symbolizes  beam  ions  that  have  not  undergone  any  collisions. 
The  energy  distribution  for  beam  ions  is  between  150  eV  and  550  eV,  with  a  mean  value  around  325  eV.  This 
suggests  a  liberal  thermal  component  that  requires  additional  investigation. 
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VI.  Results 


Computer  models  provide  an  insight  into  the  less  tangible  characteristics  of  the  plume.  Using  the  test  case 
above,  plots  were  constructed  for  several  properties  including  potential,  ion  temperature,  ion  density,  and  Debye 
length.  These  plots  are  shown  in  Figures  14-21  in  both  a  full  chamber  perspective  and  a  close  up  view  of  the 
thruster.  Potential  plots  are  shown  in  Figures  14  and  15.  As  is  predicted,  the  greatest  potential  (24  V)  is  seen 
directly  in  front  of  the  thruster,  expanding  radially  outwards  and  decreasing  with  increasing  distance  from  the 
thruster.  In  the  detail  view,  the  potential  contours  converge  along  the  thruster  centerline  a  distance  out  from  the 
thruster  face.  Ion  temperatures,  shown  in  Figures  16  and  17  follow  the  same  trend,  with  a  maximum  temperature  of 
5  eV  directly  in  front  of  the  thruster.  Ion  density,  illustrated  in  Figures  18  and  19  show  the  highest  concentration  of 
ions  (3.6  x  1017/m3)  directly  in  front  of  the  thruster,  dropping  off  rapidly  outside  the  immediate  region  of  the  plume. 
Debye  length  is  shown  in  Figures  20  and  21.  The  Debye  length  inside  the  chamber  drops  from  2.6  x  10"5  in  front  of 
the  thruster  to  8  x  10"6  at  the  back  of  the  chamber. 

VII.  Conclusion 

This  particular  test  case  has  shown  that  the  results  from  AQUILA  agree  with  experimental  data  in  magnitude  and 
trend  while  not  always  agreeing  in  distributions.  Results  from  near-field  velocity  profiles  indicate  that  the  source 
model  velocity  distributions  are  consistent  with  experimental  data,  although  results  from  the  far  field  energy 
distributions  suggest  that  the  thermal  component  may  need  additional  investigation.  The  comparisons  against 
current  density  and  energy  distribution  indicate  a  narrow  beam  divergence,  possibly  caused  by  low  collision  rates. 
Far  field  energy  profiles  also  suggest  additional  investigation  into  far  field  model  parameters,  such  as  mesh 
refinement.  Future  work  needs  to  be  completed  in  AQUILA  to  determine  which  parameters  effectively  decide  beam 
divergence,  control  velocity  distributions,  and  far  field  properties. 

VIII.  Future  Work 

The  simulation  was  completed  using  only  one,  optimal  test  case.  The  results  indicate  that  adjustments  to  the 
collision  rates,  ion  temperature  distribution,  and  far  field  properties  may  yield  more  accurate  results.  Also,  by 
changing  other  parameters  of  the  simulation  such  as  background  pressure,  mesh  quality,  and  potential  solver  inputs, 
the  results  from  the  simulation  could  change.  A  sensitivity  analysis  of  AQUILA  to  input  parameters  is  of  future 
interest.  In  addition,  more  data  from  the  BHT-HD-600  may  become  available  in  the  future.  At  that  time,  it  would 
be  pertinent  to  compare  the  data  to  results  from  this  model. 
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Figure  14:  Chamber  view  of  electrostatic  potential.  The  potential  profile  diverges  radially  outward  in  the  main 
region  of  the  beam,  decreasing  in  magnitude  from  24  V  to  0  V  at  the  back  of  the  chamber. 


Figure  15:  Thruster  view  of  electrostatic  potential.  The  potential  quickly  drops  off  from  24  V  in  front  of  the 
thruster  and  converges  along  the  centerline  of  the  thruster  downstream. 
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Figure  16:  Chamber  view  of  ion  temperature.  The  temperature  ranges  from  5  eV  in  front  of  the  thruster  to 
0.5  eV  in  the  back  of  the  chamber.  The  profile  radiates  outward  following  the  same  trend  seen  in  the 

potential. 


Figure  17:  Thruster  view  of  the  beam  ion  temperature.  The  ion  temperature  matches  the  potential  profile  by 
converging  along  the  centerline  downstream  of  the  thruster. 
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Figure  18:  Chamber  view  of  ion  density.  The  greatest  concentration  of  ion  density  is  directly  in  front  of  the 
thruster.  The  ion  density  drops  outside  the  main  region  of  the  plume  to  a  background  density  of  2  x  1016. 


Figure  19:  Thruster  view  of  ion  density.  The  ion  density  also  converges  on  the  centerline  of  the  thruster, 


similar  to  the  potential  and  ion  temperature  profiles  as  seen  in  Figures  15  and  17. 
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Figure  20:  Chamber  view  of  Debye  length.  The  Debye  length  is  1.0  x  10"5  in  the  wings  of  the  thruster  and 

drops  to  8  x  10"6  in  the  far  reaches  of  the  chamber. 


Figure  21:  Thruster  view  of  Debye  length.  The  Debye  length  in  the  front  of  the  thruster  is  2.6  x  10"5  and  3  x 

10'5  downstream,  along  the  centerline  of  the  thruster. 
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